
# Replication file for Appendix H

#install.packages("Hmisc")
require(Hmisc)
#install.packages("stargazer")
require(stargazer)

load("Brexit_Rep.RData")



eu.supp1.sub <- subset(brexit.sub2, eu.disapprove == 1)
eu.supp2.sub <- subset(brexit.sub2, eu.disapprove == 2)
eu.supp3.sub <- subset(brexit.sub2, eu.disapprove == 3)
eu.supp4.sub <- subset(brexit.sub2, eu.disapprove == 4)

dvs <- Cs(migs.take.jobs.6,migs.more.terror.6,closed.immigrants.6,hurt.standing.refugee.6,
          threaten.culture.refugee.6,overwhelm.welfare.refugee.6)

ate.1 <- vector(length(dvs), mode = "list")
names(ate.1) <- dvs

ate.2 <- vector(length(dvs), mode = "list")
names(ate.2) <- dvs

ate.3 <- vector(length(dvs), mode = "list")
names(ate.3) <- dvs

ate.4 <- vector(length(dvs), mode = "list")
names(ate.4) <- dvs


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],"~ post")
  print(modelformula)
  
  ate.1[[dvs[i]]] <- eval(substitute(lm(.modelformula, data = eu.supp1.sub, 
                                        weights = w8w6), list(.modelformula = modelformula)))
  
}


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],"~ post")
  print(modelformula)
  
  ate.2[[dvs[i]]] <- eval(substitute(lm(.modelformula, data = eu.supp2.sub, 
                                        weights = w8w6), list(.modelformula = modelformula)))
  
}

for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],"~ post")
  print(modelformula)
  
  ate.3[[dvs[i]]] <- eval(substitute(lm(.modelformula, data = eu.supp3.sub, 
                                        weights = w8w6), list(.modelformula = modelformula)))
  
}

for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],"~ post")
  print(modelformula)
  
  ate.4[[dvs[i]]] <- eval(substitute(lm(.modelformula, data = eu.supp4.sub, 
                                        weights = w8w6), list(.modelformula = modelformula)))
  
}


stargazer(ate.1[dvs[1]], ate.2[dvs[1]],ate.3[dvs[1]],ate.4[dvs[1]])
stargazer(ate.1[dvs[2]], ate.2[dvs[2]],ate.3[dvs[2]],ate.4[dvs[2]])
stargazer(ate.1[dvs[3]], ate.2[dvs[3]],ate.3[dvs[3]],ate.4[dvs[3]])
stargazer(ate.1[dvs[4]], ate.2[dvs[4]],ate.3[dvs[4]],ate.4[dvs[4]])
stargazer(ate.1[dvs[5]], ate.2[dvs[5]],ate.3[dvs[5]],ate.4[dvs[5]])
stargazer(ate.1[dvs[6]], ate.2[dvs[6]],ate.3[dvs[6]],ate.4[dvs[6]])

cis.1 <- matrix(NA,length(dvs),2)
b.1 <- matrix(NA,length(dvs),1)
cis.2 <- matrix(NA,length(dvs),2)
b.2 <- matrix(NA,length(dvs),1)
cis.3 <- matrix(NA,length(dvs),2)
b.3 <- matrix(NA,length(dvs),1)
cis.4 <- matrix(NA,length(dvs),2)
b.4 <- matrix(NA,length(dvs),1)


for (i in 1:length(dvs))
{
  b.1[i] <- ate.1[[i]]$coefficients[2]
  print(ate.1[[i]]$coefficients[2])
  ci <- confint(ate.1[[i]])
  #print(ci)
  lb <- ci[2,1]
  ub <- ci[2,2]
  cis.1[i,] <- cbind(lb,ub)
}

for (i in 1:length(dvs))
{
  b.2[i] <- ate.2[[i]]$coefficients[2]
  ci <- confint(ate.2[[i]])
  print(b.2[i])
  lb <- ci[2,1]
  ub <- ci[2,2]
  cis.2[i,] <- cbind(lb,ub)
}

for (i in 1:length(dvs))
{
  b.3[i] <- ate.3[[i]]$coefficients[2]
  print(ate.3[[i]]$coefficients[2])
  ci <- confint(ate.3[[i]])
  #print(ci)
  lb <- ci[2,1]
  ub <- ci[2,2]
  cis.3[i,] <- cbind(lb,ub)
}

for (i in 1:length(dvs))
{
  b.4[i] <- ate.4[[i]]$coefficients[2]
  print(ate.4[[i]]$coefficients[2])
  ci <- confint(ate.4[[i]])
  #print(ci)
  lb <- ci[2,1]
  ub <- ci[2,2]
  cis.4[i,] <- cbind(lb,ub)
}



y.lab <- seq(from = 1, to=length(dvs),by=1)
vnames <-c(migs.take.jobs.6 = "Migs. Take\nJobs",migs.more.terror.6 = "Migs. Bring\nTerror", 
           closed.immigrants.6 = "Not\nOpen to\nMigs.",
           hurt.standing.refugee.6 = "Refs. Don't\nImprove\nUK Image", 
           threaten.culture.refugee.6 = "Refs.\nThreaten\nCulture", 
           overwhelm.welfare.refugee.6 ="Refs.\nOverwhelm\nServices")


jitter <-  0.1
jitter2 <-  0.2
jitter3 <-  0.3

grays <- gray.colors(30, start = 0.3, end = 0.9, gamma = 2.2, alpha = NULL)


pdf("ATEs_Dis.pdf", height = 5, width=5)
par(mar =c(1.2,3.1, 4.1, 0.01), oma=c(1,5,.2,.2), cex.axis=.8, cex.main=.9)
plot(b.1, y.lab, type="n", ylab ="", xlab = "", yaxt="n", xlim = c(-0.5,0.5), ylim = c(1,6.2))
axis(2, at=y.lab, labels=vnames, las = 2, cex.axis=0.8)
#mtext("ACME",side=1,line=2,outer=F, cex = 0.8)
#segments(br.take.jobs.r$coefficients[2],x.lab, m.ci.2.lp, yr.lab, lty=1)
points(b.1, y.lab, pch= 19, cex=.8, col = "black")
segments(cis.1[,1], y.lab, cis.1[,2], y.lab, lty=1, col = "black")
segments(cis.2[,1], y.lab+jitter, cis.2[,2], y.lab+jitter, lty=1, col = grays[1])
points(b.2, y.lab+jitter, pch= 19, cex=.8, col = grays[1])
points(b.3, y.lab+jitter2, pch= 19, cex=.8, col = grays[11])
segments(cis.3[,1], y.lab+jitter2, cis.3[,2], y.lab+jitter2, lty=1, col = grays[11])
points(b.4, y.lab+jitter3, pch= 19, cex=.8, col = grays[18])
segments(cis.4[,1], y.lab+jitter3, cis.4[,2], y.lab+jitter3, lty=1, col = grays[16])
abline(v=0)
legend(0.1, 6.2, fill=c("black",grays[1],grays[11],grays[16]), 
       legend=c("Strongly\nApprove", "Approve", "Disapprove","Strongly\nDisapprove"), 
       cex =.8, border = "black", y.intersp = 1.5)
dev.off()





